Scaling of thermal conductivity of helium confined in pores 
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We have studied the thermal conductivity of confined superfluids on a bar-like geometry. We 
use the planar magnet lattice model on a lattice H x H x L with L 2> H. We have applied 
open boundary conditions on the bar sides (the confined directions of length H) and periodic along 
the long direction. We have adopted a hybrid Monte Carlo algorithm to efficiently deal with the 
critical slowing down and in order to solve the dynamical equations of motion we use a discretization 
technique which introduces errors only 0((<5t) 6 ) in the time step St. Our results demonstrate the 
validity of scaling using known values of the critical exponents and we obtained the scaling function 
of the thermal resistivity. We find that our results for the thermal resistivity scaling function are in 
very good agreement with the available experimental results for pores using the temperature scale 
and thermal resistivity scale as free fitting parameters. 
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I. INTRODUCTION 

The superfluid transition of liquid 4 He offers a unique 
opportunity for testing the finite-size scaling theory of 
static and dynamic critical phenomena. Recently, a so- 
phisticated experimental study was carried out in micro- 
gravity environment, the so-called confined helium ex- 
periment (CHeX). Lipa et al.[[l] measured the specific 
heat of helium confined in a parallel plate geometry with 
a spectacular nanokelvin resolution, thus, providing ex- 
perimental results within a few nanokelvin of T\. When 
the critical temperature is approached, the bulk corre- 
lation length £ of the fluid can become of the order of 
the confining length. CHeX approached so close to the 
lambda point that the correlation length became macro- 
scopic in size. In this case the whole fluid acts in a cor- 
related way and this changes the value of global proper- 
ties, such as the specific heat, relative to their bulk val- 
ues. In a parallel approach, Mehta and GasparinijU ||] 
have also reported earth-bound measurements on sam- 
ples with smaller plate spacing L. The size of L used 
is these measurements is smaller so that the results are 
not significantly influenced by the change in T\ between 
the top and bottom of the film because of hydrostatic 
pressure difference which exists due to the earth's gravi- 
tational field. 

The finite-size scaling (FSS) theory [Q and the 
renormalization-g roup theory (RGT)§ were expected to 
describe the behavior of the system at temperature close 
to T\. A testable implication of this theory is that very 
close to the lambda point, in a confined system with a 
confining length of size H, a dimensionless quantity or 
the ratio of two quantities having the same dimensions, 
is only a function of the ratio oi^/H. Therefore the val- 
ues of a given observable 0(t, H), for various values of H 
and of the reduced temperature t = \1 —T/T\\, divided 
by its bulk value of 0(t, H = oo) should be a dimension- 
less scaling function f(x), where x = £(t)/H. The re- 
sults of CHeX were in remarkable agreement with predic- 
tions which were available prior to the experiment based 



on scaling functions obtained from renormalization-group 
theoryQ and those obtained by combining FSS and the 
results obtained from large-scale simulations||]. 

A second equally important step toward understand- 
ing the FSS theory is to study dynamical and transport 
properties near a critical point. A well-suited candidate 
problem for this study is the thermal conductivity A of 
4 He near T\. When T\ is approached from above, the 
thermal conductivity of the fluid diverges 111]]. The 
precise behavior of bulk A as a function of t was stud- 
ied in great detail both experimentally |l2|, 13, |l4j and 
theoretically fl5||. 

There are several recent theoretical studies of dynam- 
ical critical phenomena and dynamical scaling. Koch, 
Dohm, and Stauffer|16 presented field-theoretical and 
numerical studies of the validity of dynamic finite-size 
scaling for relaxational dynamics in cubic geometry with 
periodic boundary conditions above and below T c . Quan- 
titative agreement between theory and Monte Carlo data 
was obtained by them. Koch and DohmJlTj have pro- 
vided a prediction for the dynamic finite-size scaling func- 
tion for the effective diffusion constant of model C of Ho- 
henberg and Halperin|l^]. Bhattacharjee|l9t derived an 
approximate form of the scaling function for the thermal 
conductivity using a decoupled-mode approximation and 
model E. Krech and Landau pf| calculated the transport 
coefficient of the out-of-plane magnetization component 
at the critical point, which is related to the thermal con- 
ductivity of liquid He using Monte Carlo spin dynamics 
simulations of the XY model in three dimensions on a 
simple cubic lattice with periodic boundary conditions. 
They determined the critical exponent characterizing the 
thermal conductivity. 

Accurate experimental studies have been carried out 
not only for dynamic bulk phenomena with improved 
resolution but also dynamic properties in confined ge- 
ometries deeply in the critical region pl| p2[ . Rather re- 
cently, Kahn and Aiders |2j| measured the thermal con- 
ductivity of liquid 4 He confined in a glass capillary array 
of thickness 3 mm with holes 2/im in diameter. Their 
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results show that long cylindrical samples have a transi- 
tion from three-dimensional to one-dimensional behavior 
and there is no phase transition in the one-dimensional 
system. However, as measurements over a wide range of 
hole-diameter are required in order to test the finite-size 
scaling theory for transport properties, further experi- 
mental studies are planned p4[ in order to reveal dynam- 
ical exponents near the critical point and to study the 
finite-size scaling behavior of the thermal conductivity in 
such confining geometries. To avoid the limitations im- 
posed by the earth's gravity, this experimental effort 1 24 
will be carried out under microgravity conditions on the 
Low Temperature Microgravity Facility on International 
Space Station. 

In this paper we wish to study the thermal conduc- 
tivity A of confined helium and to calculate the scaling 
function associated with A for a fixed geometry. Since 
there are already experimental results [|23| for the scaling 
function of A for the pore-like geometry, in this paper 
we will focus our attention to this geometry because we 
hope to compare with the experiment. We will examine 
the FSS theory for the thermal conductivity of helium 
confined in a bar-like geometry i.e., on an H 2 x L lat- 
tice with L >> H. This confining geometry is similar 
to that of Kahn and Ahlers[^3| because two of the di- 
mensions of a pore used in their experimental studies are 
confining as is the case of the bar-like geometry. We will 
consider the limit in which our results arc independent of 
the bar length L. We will apply periodic boundary condi- 
tions (BC) in the L direction because these BC approach 
the bulk limit faster. In the other two directions which 
are kept finite we will apply open boundary conditions. 
We will use the dynamics of planar-magnet model and 
Monte Carlo simulation to study X(t,H). We find that 
A(t, H)H- R I V plotted as a function of x = tH » fall on 
the same curve for a wide range of values of H and t, us- 
ing the known values of v and 7r. This demonstrates that 
finite-size is also valid for dynamical critical properties. 
In addition we obtain the scaling function which fits very 
well the experimental data of Kahn and Ahlers[^3) using 
the scale of temperature and the thermal conductivity 
scale as free parameters. 



II. THE METHOD 

We will first briefly describe the model and the numer- 
ical method used and show how the thermal conductivity 
is computed in our model. To describe the dynamics of a 
superfluid, we will use the planar magnet model which is 
classified as model F (or E in the absence of an external 
field) by Hohenberg and Halperin|ll|. Matsubara and 
Matsuda[^5| has proposed model F to explain the prop- 
erties of liquid 4 He. The problem of hard core bosons can 
be described by a lattice gas model which can be mapped 
to the quantum antiferromagnet in which the superfluid 
order parameter corresponds to S x — iS y while the den- 
sity of the boson system corresponds to 1/2 — S z . In 



order to study equilibrium critical properties of a super- 
fluid one uses the XY model || because the planar 
magnet model and the XY model belong to the same 
universality class (2^. For critical dynamics of a super- 
fluid, however, one needs to use the full planar magnet 
model in which the role of the third component of the 
pseudospin is crucial. [fl8| 

In the pseudospin notation, the planar magnet model 
takes the following form: 



H=-J^(SfS* + S?S%), 



(1) 



where the summation is over all nearest neighbors, Si — 
(Sf , Sf , Sf), and J sets the energy scale. In the usual 
XY model the two-component pseudo-magnetization 
corresponds to the superfluid order parameter. In the 
planar magnet model, the third component corresponds 
to the particle density and it is necessary in order to 
study the dynamics. 

In our calculations, we use a bar-like geometry, i.e. a 
HxHxL lattice with L H . This geometry is chosen in 
order to mimic the pore geometry used in experimental 
studies. In our calculations, open boundary conditions 
are imposed in the H direction, and in the L direction 
we applied periodic boundary conditions. 

We use a hybrid Monte Carlo procedure [^oj which 
consists of a combination of steps using the Metropolis 
update, the Cluster update^, and the over-relaxation 
algorithm ]29{. Using this hybrid algorithm, we generate 
approximately 3,000-10,000 uncorrelated configurations 
from the equilibrium canonical ensemble at a given tem- 
perature. Each configuration is evolved using the equa- 
tions of motion for the planar magnet model which are 
given as follows (TH, [20| 



dt 



dH 



x Si 



(2) 



Starting from a particular initial spin configuration, we 
perform numerical integration of these equations of mo- 
tion. Following Ref. ||3C|l we use a recently developed de- 
composition method ||31| where the integration is carried 
out to a maximum time t max (typically of the order of 
£maa;=400) with a time step <5i=0.05. We made sure that 
this way we determined the real-time history of every 
configuration within a sufficiently long interval of time 
(0< t < t max ). Finally, we compute the thermal aver- 
age of a time-dependent observable (such as the thermal 
current-current correlation function) by averaging over 
all the values of the observable obtained by evolving all 
the independent initial equilibrium configurations gener- 
ated via the hybrid Monte Carlo procedure. 

Compared to calculating static critical properties, the 
computation of dynamical properties is far more CPU 
time intensive and demands large computational re- 
sources. The computations described here were carried 
out on a dedicated massively parallel cluster of 64 CPUs 
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which was designed by our group to achieve high perfor- 
mance to cost ratio. 

We computed the thermal conductivity on H x H x L 
lattices, where #=6,8,10,12,14,20 and L = 5H. The 
thermal conductivity A of liquid 4 i/e at a given temper- 
ature T can be calculated using the dynamic current- 
current correlation function p0|: 

A=tt^ / dtJ2<3o(0)3?(t)>, (3) 

kbTxzz tt Jo 

where the out-of-plane static susceptibility 

Xzz =< M z2 > /(k B TL 3 ) (4) 

is needed for normalization. The ^-component j* of the 
current density ji associated with the lattice point i is 
defined by 

3i = J{SiSi+e z ~ ^i^i+e-^i (5) 

where the notation i + e z denotes the nearest neighbor of 
the lattice site i in the z lattice direction. 

Now, we would like to examine the finite-size scal- 
ing hypothesis for the thermal resistivity R(t, H) = 
l/X(t,H), and to compare our results with the existing 
experimental results |23j. The dependence upon t of the 
bulk thermal resistivity can be described by the power 
law 

R(t) = Rot 71 , (6) 

where ir is a dynamic critical exponent. Using Eq. |^, the 
finite-size scaling expression for the thermal resistivity 
R(t, H) is given by 

R{t,H)H w ' v = f{tH^ u ), (7) 

where the function f(x) is universal and ^=0.6705 is the 
critical exponent of the correlation length |32|. 

III. RESULTS 

In this section, we calculate the thermal resistivity, we 
examine its scaling behavior with respect to H and then 
we compare the scaling function with the experimental 
results. To calculate these observables with small statis- 
tical errors even with our utilization of the most recent 
numerical advances and with using the 64-node dedicated 
cluster, it requires signinificant amount of time of high- 
throughput computation. 

Fig. [I] shows some of our results for the thermal resis- 
tivity R{T, H) as a function of temperature T for var- 
ious lattice sizes with open boundary conditions in the 
H direction. Our results for R(T, H) for several values 
of H and T are given in Table § and 0. We wish to 
make L large enough so that finite-size effects due to L 
are smaller than our statistical errors. We have found 
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FIG. 1: Thermal resistivity R(T,H) versus temperatures for 
bar-like lattices with sizes that correspond to H = 6, 8, 12, 20 
and L = 5H. The bulk T A = 1.5518 is also shown. 

10 i . 1 . 1 . 1 . 




-5 5 10 15 



FIG. 2: The universal function f(x) obtained for bar-like ge- 
ometry. The solid line represents the available experimental 
results for pore-like geometry. In the experimental results the 
resistivity scale and the temperature scale are used as free 
parameters. 



that taking L « 5H and applying periodic boundary 
conditions along the direction of L introduces insignif- 
icant finite size effects due to the finite size of L for the 
temperature range studied here. Since we wish to re- 
main in the 3D critical region, it appears that keeping 
L = 5H introduces small finite-size effects due to L in 
this region. Lowering the temperature further, when the 
correlation length becomes comparable to L, the value of 
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FIG. 3: The universal function f(x) obtained for bar-like ge- 
ometry using the theoretical value of tt = 0.335. 



R is influenced by the finite size of L. In this region the 
system exits the 3D critical region and it behaves as a 
one-dimensional system. 

Notice in Fig. |l|that the thermal resistivity feels strong 
finite-size effects due to the bar thickness H . The ar- 
row shows the bulk transition temperature T\ — 1.5518 
obtained from Monte Carlo simulation using the planar 
magnet model[^. In bulk helium R(t) approaches zero 
as the bulk transition temperature T\ is approached from 
above. 



T/J 


H = 6 


H = 8 


H = 10 


# = 12 


H = 14 


1.40 


0.350(34) 


0.178(11) 


0.053(4) 


0.033(4) 




1.45 


0.535(35) 


0.286(17) 


0.143(10) 


0.068(6) 


0.050(4) 


1.50 


0.662(61) 


0.470(36) 


0.353(38) 


0.247(31) 


0.182(21) 


1.5518 


0.843(62) 


0.670(56) 


0.614(53) 


0.501(56) 


0.452(59) 


1.60 


0.923(84) 


0.821(36) 


0.688(67) 


0.652(55) 


0.706(42) 


1.65 


1.028(81) 


0.951(64) 


0.854(42) 






1.70 




1.114(69) 


0.901(63) 


0.988 (102) 


0.984(105) 


1.80 


1.213(97) 


1.216(86) 




1.125(94) 


1.081(166) 



TABLE I: Calculated results for the thermal resistivity for 
lattices H x H x L with L m 5H and H = 6.8, 10, 12, 14. 
The number in parenthesis gives the error in the last decimal 
places. 



We wish to avoid using any adjustable parameters to 
obtain scaling of our results. Thus, we need to examine 
if our results obey scaling using the known values of the 
critical exponents v and tt. The value of v is accurately 
known from theoretical and experimental studies of static 
critical properties and we shall use the value v — 0.6705 



T/J 


H = 20 


1.50 


0.053(3) 


1.52 


0.158(13) 


1.54 


0.294(27) 


1.56 


0.408(48) 


1.58 


0.572(79) 


1.60 


0.567(62) 


1.65 


0.733(76) 


1.75 


1.086(130) 



TABLE II: Calculated results for the thermal resistivity for 
an 20 x 20 x 100 size lattice. 



as determined by Goldner and Aiders j32|. There is less 
agreement between theory and experiment on the actual 
value of the dynamical critical exponent tt. Ahlers|33jl 
used a power law fit to the data of Tarn and AhlersjlJ] 
for their "Cell F" and he found the value tt = 0.4397. 
However, the dynamic scaling theory had predicted 
a divergence in A with a critical exponent given by tt — 
i//2» 0.335. 

Fig. H shows a scaling plot of the thermal resistivity 
scaling function f(x) — R(t, H)t n / V versus the scaled 
reduced temperature parameter x = tH~ , where the re- 
duced temperature is taken relative to the bulk transition 
temperature T\. Our Monte Carlo data collapse onto a 
universal curve using the value of tt w 0.44 determined 
by Aiders j33|. In Fig. || we compare our universal func- 
tion fix) with the experimental data obtained by Kahn 
and Ahlers[^3| represented by a solid line. In order to do 
this, we used two multiplicative constants as free fitting 
parameters, one multiplying the scale of x axis and an- 
other the scale of y. The agreement between Monte Carlo 
simulation and experiment is quite satisfactory. In the 
past it has been demonstrated |2(| [55| that the boundary 
conditions play a significant role in defining the universal 
function f(x). We believe that if we use more realistic 
boundary conditions, such as Dirichlet boundary condi- 
tions, along the 7?-direction we can reduce the number 
of fitting parameters to only one. 

Using the theoretical value of n, the results of our sim- 
ulation also collapse on a different scaling function given 
in Fig. |3| . However, if we attempt to fit the scaling curve 
with the experimental resistivity of Kahn and Ahlers we 
obtained a lower quality fit than that of Fig. |2|. 

In summary we have calculated the thermal resistiv- 
ity R(t, H) of liquid 4 if e in a pore-like geometry (on a 
H x H x L lattice) applying open boundary conditions in 
the H direction. We have been able to demonstrate the 
validity of finite-size scaling theory and we obtained the 
thermal resistivity scaling function f[x) using known val- 
ues for the critical exponents and no adjustable param- 
eters. In addition, the scaling function f[x) for R(t, H) 
agrees rather well with experimental data using the tem- 
perature scale and thermal resistivity scale as free pa- 
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rameters. 
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